Computation of the response functions of spiral waves in active media 
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Rotating spiral waves are a form of self-organization observed in spatially extended systems of physical, 
chemical, and biological nature. A small perturbation causes gradual change in spatial location of spiral's 
rotation center and frequency, i.e. drift. The response functions (RFs) of a spiral wave are the eigenfunctions 
of the adjoint linearized operator corresponding to the critical eigenvalues A = 0, ±iuo. The RFs describe the 
spiral's sensitivity to small perturbations in the way that a spiral is insensitive to small perturbations where 
its RFs are close to zero. The velocity of a spiral's drift is proportional to the convolution of RFs with the 
perturbation. Here we develop a regular and generic method of computing the RFs of stationary rotating spirals 
in reaction-diffusion equations. We demonstrate the method on the FitzHugh-Nagumo system and also show 
convergence of the method with respect to the computational parameters, i.e. discretization steps and size of the 
medium. The obtained RFs are localized at the spiral's core. 

PACS numbers: 02.70.-c, 05.10.-a, 82.40.Bj,82.40.Ck, 87.10.-e 



I. INTRODUCTION 

Autowave vortices, or spiral waves in two-dimensions (2D), 
are types of self-organization observed in dissipative media of 
physical [1-4], chemical [5-7], and biological nature [8-13], 
where wave propagation is supported by a source of energy 
stored in the medium. The common feature of all these phe- 
nomena is that they can be mathematically described, with 
various degrees of accuracy, by reaction-diffusion partial dif- 
ferential equations, 

d t u = f (u) + DV 2 u, u,fel', D G R txt , i>2, (1) 

where u(r, t) = (ui, . . -ui) T is a column-vector of the 
reagent concentrations, f(u) = (fi, . . . fe) T is a column- 
vector of the reaction rates, D is the matrix of diffusion coef- 
ficients, and f £ R 2 is the vector of coordinates on the plane. 

The existence of vortices is not due to singularities in the 
medium but is determined only by development from initial 
conditions. A rigidly rotating spiral wave solution to the sys- 
tem ( 1 ) has the form 

U = \J(p(f-R),$(f-R)+ut-<f>), (2) 
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where p(r — R), $(r — R) are polar coordinates centered at 
R, vector R = (X, Y) T defines the center of rotation, and $ 
is the initial rotation phase. For a steady, i.e. rigidly rotating, 
spiral R and <I> are constants. The system of reference co- 
rotating with the spiral's initial phase and angular velocity to 
around the spiral's center of rotation is called the system of 
reference of the spiral. In this system of reference, R = 0, 
$ = 0, and the polar angle is given by 9 = -d + ut. In this 
frame the spiral wave solution U(p, 6) does not depend on 
time and satisfies the equation 

f (U) - cjU 9 + DV 2 U = 0. (3) 

In this equation, the unknowns are the field U(p, 9) and the 
scalar w. 

A slightly perturbed steady spiral wave solution 

V(p,9,t)=XJ(p,6)+e S (p,6,t), gel', 0<e«l, 

substituted in (1), at leading order in e, yields the evolution 
equation for the perturbation g, 

a tg = a u f(u)g-^g + Dv 2 g. 

Thus, the linear stability spectrum of a steady spiral 

CV = XV (4) 
is defined by the linearized operator 

£ = DV 2 - Lud e + <9 u f (U). (5) 
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The operator L has critical (Re (A) = 0) eigenvalues 

A„ = inuj, n = 0, ±1, (6) 

which correspond to eigenfunctions related to equivariance of 
(1) with respect to translations and rotations, i.e. "Goldstone 
modes" (GMs) [14-17] 

V<°> = -doUQ>,6), 

V {±1) = -^e* (dpTip-^VfaO). (7) 

The stability spectra of steady spiral waves was originally ob- 
tained numerically by Barkley [16]. Subsequently the spec- 
trum was analysed for infinite and large bounded domains by 
Sandstede and Scheel [18-20] with follow-on numerical in- 
vestigations by Wheeler and Barkley [21] confirming the large 
domain behavior of the stability spectrum. 
In a slightly perturbed problem 

d t u = f(u) +DV 2 u + eh, hel', < e < 1, (8) 

where eh(u, r, t) is some small perturbation, spiral waves 
may drift, /.e.change rotational phase and/or center location. 
Then, the center of rotation and the initial phase are no longer 
constants but become functions of time, R = R(t) and 
$ = $(*). 

In linear approximation, assuming that 

k, $ = 0(e), 
the drifting spiral wave solution can be represented as 

U = \J(p{r-R(t)), d(?-R(t))+wt-$(t))+eg(?, t), (9) 

where eg(r, t) is a small perturbation of the steady spiral 
wave solution U. 

Then, the solution perturbation g in the laboratory frame of 
reference will satisfy the linearized system 

(d t - Dv 2 -a u f(u)) g 

= h(u,F,t) - -(R- V + <f>d e )U. (10) 

e 

The solvalability condition for equation (10) for g , i.e. Fred- 
holm alternative, re-written in the spiral frame of reference, 
requires that the free term must be orthogonal to the kernel 
of the adjoint operator to C defined in (5). This leads to the 
following system of equations for the drift velocities 

$ = eF (R,t), R = eFx(R,t). (11) 

Thus, the drift velocities $ and R are determined by the 
"forces" F Q and Fx = (Re (Fx) , Im (Fx)) T which, after slid- 
ing averaging (more specifically, central moving average) over 
the spiral wave rotation period, can be expressed [15] as 

F n (R,t) = e ira * j> 

t — -K I UJ 

x (V ( ' l) (p(f- R), d(f- fl)+ur-$) , h(f, t)\ , 
n = 0,±l. (12) 



(of course, F_i = Fx). Here (■ , •) stands for the scalar prod- 
uct in functional space, 



(w , v) = 




r- 



The kernels W^") of convolution-type integrals in (12) are 
the spiral wave's response functions (RFs), i.e., the critical 
eigenfunctions 

£+W ( "> = n n W (n \ (13) 

where 

fi n = —iivn, n = 0, ±1, (14) 
of the adjoint linearized operator: 

£+ =T>V 2 +ujdo + (d u i(U)) T , (15) 
chosen to be biorthogonal 

(w^,VW)=^, fcj (16) 

to the Goldstone modes (7). Note that the RFs do not depend 
on time, i.e. are functions of the coordinates only, in the co- 
rotating system of reference. 

The asymptotic theory just outlined reduces the description 
of the smooth dynamics of spiral waves from the system of 
nonlinear partial differential equations (1) to the system of or- 
dinary differential equations (11), describing the movement 
of the core of the spiral and the shift of its angular velocity. 
Several qualitative results in the asymptotic theory of spiral 
and scroll dynamics have been obtained without the use of re- 
sponse functions, e.g. [15, 17, 22-30]. However, an explicit 
knowledge of RFs makes possible a quantitative description, 
which obviously can be much more efficient for the under- 
standing and control of spiral wave dynamics in numerous ap- 
plications, e.g. control of re-entry in the heart. 

The asymptotic properties of the RFs at large distances 
are crucial for convergence of the convolution integrals in 
(12). An early version of the asymptotic theory, developed 
by Keener [31] for scroll wave dynamics, considered the RFs 
asymptotically periodic in the limit p — > oo, in much the same 
way as spiral waves are, thus requiring an artifical cut-off pro- 
cedure to tackle the divergence of the integrals in (12) follow- 
ing from such an asumption. 

Based on observations and empirical data of spiral wave 
dynamics, Biktashev [14, 32] conjectured that the response 
functions quickly decay at large p, i.e. are effectively local- 
ized. This conjecture implies that the integrals in (12) con- 
verge and no cut-off procedure is required. 

To prove existence of the localized responce functions, Bik- 
tasheva et al. [33] explicitly computed them in the complex 
Ginzburg-Landau equation (CGLE) for a particular set of pa- 
rameters. Those computations exploited an additional sym- 
metry present in the CGLE, which permitted the reduction of 
the 2D problem to the computation of ID components. The 
computations were verified by numerical convergence of the 
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method with respect to the space discretisation and the size 
of the medium. Following this work, the computed RFs were 
successfully used for quantitative prediction of the spiral's res- 
onant drift and drift due to media inhomogeneity [34, 35]. By 
explicitly computing the RFs in the CGLE for a broad range 
of the model's parameters, Biktasheva and Biktashev [36, 37] 
showed that the RFs are localized for stable spiral wave solu- 
tions and qualitatively change at crossing the charachteristic 
lines in the model parameter plane. 

Recently, there has been a significant theoretical progress 
in mathematical treatment of the localization of the response 
functions. Sandstede and Scheel [38, Corollary 4.6] analyti- 
cally proved such localization for one-dimensional wave dis- 
locations, which may be considered as analogues of a spiral 
wave in one spatial dimension. Hopefully this can be extended 
to two spatial dimensions, i.e. to spiral waves. 

For cardiac applications, dynamics of spiral waves in ex- 
citable media is more important than in oscillatory media such 
as the CGLE, as most cardiac tissues are excitable. These 
models do not allow reduction to ID, making quantitatively 
accurate computation of the response functions more chal- 
lenging. So far, the response functions have been computed 
in the Barkley [39, 40] and FitzHugh-Nagumo [41] models 
of excitable media. For the chosen sets of model parameters, 
the computed RFs appeared effectively localized in the vicin- 
ity of the spiral wave core. Hamm [39] and Biktasheva et 
al. [41] calculated RFs on Cartesian grids, but the accuracy 
was not sufficient for quantitative prediction of drift. Hakim 
and Henry [40] took the advantage of a polar grid and Barkley 
model to compute the spiral wave solution with an accuracy 
of 10~ 8 and RFs with accuracy 10~ 6 (both in the sense of 
?2-norm of the residue of the discretized equations) leading to 
quantitative prediction of drift velocities with about 4% accu- 
racy. 

Encouraging as these results are, there is a need for a more 
computationally efficient, accurate and robust method to com- 
pute the response functions of spiral waves in a variety of ex- 
citable media with required accuracy. The aim of this paper 
is to present a method which is superior to previous methods 
used to compute response functions and to demonstrate that it 
works for stationary rotating spirals in FitzHugh-Nagumo sys- 
tem. We also demonstrate convergence of the method with re- 
spect to the computational parameters, i.e. discretization steps 
and size of the medium, and show that the method is vastly 
more efficient than the methods used before [40, 41]. 



II. METHODS 
A. Computations 

To compute the response functions, we use methods similar 
to those described in [16, 21]. 

The nonlinear problem (3) is considered on a disk p < 
Pmax, with homogeneous Neumann boundary conditions, 
dpU(p max , 9) = 0. The fields are discretized on a regular 
polar grid (pj,9k) = (jAp,kA9) where < j < N p and 
< k < Ng plus the center point p = 0. Hence there are 



NpNg + 1 grid points and correspondingly N = £(N p Ng + 1) 
unknowns and the same number of equations in the discretiza- 
tion of (3). For the inner points j < N p , the p-derivatives 
are calculated via second-order central differences. The 9- 
derivatives are calculated using Fornberg's weights . f sub- 
routine [42] which uses all Ng values so, in theory, provides 
an approximation of ^-derivatives of the order of Ng. The dis- 
cretization of the Laplacian at the center point is via the differ- 
ence between the average around the innermost circle p = Ap 
and the center point, and the approximation at j = N p takes 
into account the boundary conditions at p = p max . 

The discretized nonlinear steady-state spiral problem (3) 
is solved by Newton's method, starting from initial approx- 
imations obtained by interpolation of results of simulations of 
the time-dependent problem (1) using EZSPIRAL. The New- 
ton iterations involve inversion of the linearized matrix which 
has a banded structure with the bandwidth 1 + 2£Ng. This is 
achieved by the appropriate ordering of the unknowns of the 
discretized problem within the TV-dimensional vector of un- 
knowns, so that the index enumerating components of reagent 
vectors from R £ varied fastest, followed by the index enumer- 
ating angular grid points k A9, followed by the index enumer- 
ating the radial grid points jAp. 

The thus posed discretized nonlinear problem inherits the 
symmetry of (3) with respect to rotations. To select a 
unique solution out of a family of solutions generated by 
this symmetry, we impose a "pinning condition" of the form 
Ue r (j*Ap, k*A9) = m», where £„, it* and j* may be selected 
arbitrarily and fc* is chosen as the #-grid point in the p = j*Ap 
circle that gives the I* -component value closest to u„ in the 
initial approximation. Since Ui, (j* Ap, fc* A9) is fixed, it is 
no longer an unknown, and its place in the R w -vector of un- 
knowns is taken by lo, also to be found from (3). In this way, 
the balance of the unknowns and equations is preserved. As 
oj is present in all equations, the corresponding non-zero col- 
umn of the linearization matrix destroys the handedness of the 
matrix. This obstacle is overcome by employing the Sherman- 
Morrison formula [43] to find solutions of the corresponding 
linear systems using only banded matrices. Newton iterations 
are performed until the residual in solution of the discretized 
version of equation (3) becomes sufficiently small. 

The linearized problems (4) and (13) are considered in the 
same domain with similar boundary conditions. The critical 
eigenvalues and eigenvectors of the discretized operators C 
and C + are computed with the help of a complex shift and 
Cayley transform. 

For a matrix L, be it discretization of C or £ + , the complex 
shift is defined as 

A = L + ikI 

and the subsequent Cayley transform as 

B = (ei + A)- 1 ( ?? I + A) (17) 

where k, £ and ?/ are real parameters and I is the identity ma- 
trix. If A, a and (3 are eigenvalues of L, A and B, respectively, 
this implies 

a = A + ik, p = . 

4 + a 
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The selected eigenvalues and eigenvectors of the thus con- 
structed matrices B are then found by the Arnoldi method, 
using ARPACK [44]. 

We have used £ = 0, r\ = 1 and k = 0, ^fui when seeking, 
respectively, V^°' ±:L ) and W^'^ 1 ), where uj is the solution 
of the corresponding nonlinear problem previously obtained. 
With this choice of £, rj and k, the numerical eigenvalues A 
and fi closest to the theoretical critical eigenvalues (6) and (14) 
correspondingly, generate the largest \ j3\. Hence, the Arnoldi 
method in each case is required to obtain the eigenvalue with 
the largest absolute value. 

To normalize the eigenvectors, we use the "analytical" 
Goldstone modes Vw > obtained by numerical differentiation 
of the numerical spiral wave solution U, namely, 

V(°) = -detiQ>,6), 
V {±1) = -le^ e (d pT ip- 1 de)lJ(p,e), 

where differentiation has been implemented using the same 
discretization schemes as used in calculations. 

First, the response functions \V i - k > computed by ARPACK 
are normalized with respect to the "analytical" Goldstone 
modes V' fe ) so that 



w« , V (fc) ) = 1, fc = 0,±l, 



where numerical integration involved in (• , •} has been carried 
out using the trapezoidal rule. 

Then, the "numerical" Goldstone modes V( fc ) computed by 
ARPACK are normalized with respect to the normalized re- 
sponse functions so that 



W (fc) ,V w =1, fc = 0,±l. 



Thus, we finally obtain 

• a numerical solution for the spiral wave problem (3) to- 
gether with the angular velocity u), 

• "analytical" Goldstone modes V' 4 ', 

• normalized "numerical" Goldstone modes V'*', and 

• normalized response functions W( fe ) . 

B. Analysis 

To validate the computed response functions, we have to 
demonstrate convergence of the solution with respect to the 
numerical approximation parameters such as the size of the 
medium p max , and the discretization steps Ap and A9. 

First of all, we have to demonstrate convergence of the 
computed eigenvalues of A„ and fi n to their theoretical values 
(6) and (14), taking for ui its numerical approximation uj found 
by numerical solving the discretized problem (3). Since the 
"theoretical" value for uj is not available, we can only check 
convergence of uj to some limit. 



The accuracy of the "numerical" Goldstone modes is quan- 
tified by the distance between the "numerical" and "analyti- 
cal" Goldstone modes, in norm 



1/2 



as well as Cn norm 



V^(r)- V^(r) 



dV 



V; = max 



r£S 



V^(r) - V^(f) 



over a disk S of half the radius of the computational domain: 

S = {f: |r| </w/2}. 

The smaller disk is used to exclude the effects of boundary 
conditions. The issue is that the exact GM V do not satisfy 
Neumann boundary conditions whereas V do, hence there is 
an inevitable deviation between them near p = p max , which 
is an artefact of restricting our problem to a finite domain, and 
is not indicative of the accuracy of the computed W, which 
are expected to be exponentially small near p = p max . 

The accuracy of the computed response functions W could 
be tested directly in the same way as the accuracy of the com- 
puted uj, i.e. by the numerical convergence to some limit. This 
is however, difficult to implement for the numerical solutions 
obtained on different grids. Nevertheless, we are able to ex- 
amine the convergence in Ap where coarser grids are subgrids 
of the finer grids by restricting the fine-grid solutions to the 
coarse grid, without the need for any interpolation. Specifi- 
cally, we calculate 



1/2 



W«,(r)-WXMr) 



U) 



Ap 



Ap, 



and 



£' 



max 

feB 



Ap 



(J) 



Ap, l 



over the whole computational domain 



B = {r:\r\< p max }, 

where W^],(r) are the numerical response functions calcu- 
lated at the radius step Ap which is an integer multiple of the 
minimal radius step Ap*, and the finest numerical response 

functions (r) have been restricted to the coarser grid 

of w2p(r) of the solution to which they are compared, so 
the numerical integration is done over the coarser grid. Note 
that in the series with varying p max and fixed AO and Ap, the 
coarser grids are also subgrids of the finer grids, but as the 
pinning point is defined via p max , solutions at different p max 
are again not directly comparable to each other so this series 
is not used in this comparison. 
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We also assess accuracy indirectly via the bi-orthogonality 
between the response functions and the Goldstone modes re- 
quired by (16). Specifically, we examine the orthogonality of 
the RFs to the "analytical" GMs, quantified by 



O a 



E E 

j=0,±l fc=0,±l 



(18) 



and orthogonality of the RFs to the "numerical" GMs quanti- 
fied by 

On= E E |(w«,VW}-^, fe | 2 . 
j=0,±l fe=0,±l 

Note, that by construction the diagonal elements of both the 
"numerical" and "analytical" bi-orthogonality matrices here 
are all equal to 1 up to round-off errors. 

The measures O a and O n require some discussion. The bi- 
orthogonality should be exact for exact RFs and GMs. How- 
ever, what we calculate are approximations of these functions, 
subject to discretization in p and 9 and restriction to a finite 
domain p < p max . The bi-orthogonality of numerical solu- 
tions is therefore not exact and its deviation from the ideal is 
an indication of the accuracy of calculation, and its conver- 
gence in Ap, A9 and p max is an indication, albeit indirect, of 
the accuracy of the solutions. 

In more detail, if the the matrices representing discretiza- 
tion of C and C + were transposes of one another, then their 
eigenvectors corresponding to different eigenvalues would be 
exactly orthogonal in 1%, and so a measure of their orthog- 
onality would not depend on the spatial discretization but 
only on the accuracy of the calculation of the eigenvectors by 
ARPACK. However, C and C + are conjugate with respect to 
the scalar product which is approximated by a discrete inner 
product with a weight, hence the matrices of C and C + are 
not transposed. Moreover, because of the approximation used 
for these operators {e.g. high-order approximation in A9 vs 
second-order approximation in Ap), the corresponding matri- 
ces are not adjoint of each other with respect to the weighted 
1-2 either. So, O n provides a measure of the consistency of 
these matrix representations together with the accuracy with 
which the eigenvectors are computed with ARPACK. 

Moreover, apart from the question of accuracy of finding 
the eigenvectors of the discretized operators and accuracy of 
finding the eigenfunctions of the original continuous opera- 
tors, there remains a question of whether the found eigenvec- 
tors and eigenfunctions are the ones that we need, that corre- 
spond to and ±iu, rather than eigenfunctions corresponding 
to eigenvalues which happened to be close to and ±ioj [49]. 
For the GMs, the answer to this question is ensured by check- 
ing the distance Vj ; however, this answer is not absolute as the 
comparison is made only over part of the disk, for reasons dis- 
cussed above. We note, however, that the C + eigenfunctions 
corresponding to the eigenvalues close to but different from 
0, ±iu>, are orthogonal to the GMs and for them O a would be 
not small [50]. Since O a is defined in terms of scalar products 
with the mode determined directly from the underlying spiral 
wave, its smallness provides the additional assurance that the 
adjoint eigenfunctions are indeed the RFs that we are after, 
not just some adjoint eigenfunctions. 



in. RESULTS 
A. General 

We have tested our method for computing the response 
functions in the case of the FitzHugh-Nagumo model, t = 2, 

fi = e^Oi - uf/3 - u 2 ), 
h = £ ( u i - au 2 + b), 

, with parameters a = 0.5, b = 0.68, e = 0.3. 



D 



For pinning, we have used =2,u* = 0.1 and = N p /2. 
Newton iterations have been performed until the Euclidean 
(I2) norm of the residual in the discretized nonlinear equa- 
tion falls below 10~ 8 . For comparison, we have also run 
cases, discussed later in fig. 5, in which iterations continue 
until the norm of the residual no longer decreases (typically 
such norms were below 10~ 9 down to 10 -13 ). The tolerance 
in ARPACK's routines znaupd and zneupd has been set 
to the default "machine epsilon". For the Krylov subspace 
dimensionality we have tried 3 and 10, with no perceptible 
difference in either the numerical results. 

Before discussing the performance of our numerical tech- 
niques, we briefly present typical solutions. Figures 1 and 2 
illustrate the spiral wave solution and the GMs and RFs for 
Pmax = 25, N p = 1280 and N 9 = 64. This solution is 
taken as the best achievable given memory restrictions (4Gb 
of real memory). The angular velocity for it was found to be 
lj w 0.5819341748776017. For the GMs and RFs, we show 
the n — and n = 1 modes only, since the calculated n = — 1 
modes are almost exactly the complex conjugates of the n = 1 
modes, which of course they should be. 

One can see that the GMs V are indeed proportional to cor- 
responding derivatives of the spiral wave solution U, and that 
the RFs W are localized in a small region of the spiral tip and 
are indistinguishable from zero outside that region. 

The character of the RFs' decay with distance is illustrated 
in more detail in fig. 3. We plot the angle-averaged values of 
the solutions, defined as 

for X = U, V and W. Note the difference in the behavior 
of (U)\ n ^ and (V) on one hand and (W)j on the other 
hand. In the semilogarithmic (linear for horizontal axis, log- 
arithmic for vertical axis) coordinates of fig. 3(c) the graphs 
of {W}^ (p) are straight for a large range of p, not too close 
to or /9 max = 25, and for several decades of magnitude of 
(W)^ . This shows clearly the expected exponential localiza- 
tion of the RFs. For comparison, we also show convergence 
of uj — cD(p max ) in a. disk as a function of the disk radius 
p max . Theory [36, 45^17] predicts that the (W)^ (p) and 
Aw(p max ) = c2>(p m ax) ~ oj(oo) dependencies should both be 
decaying exponentials with the same characteristic exponent; 
this agrees well with the numerical results shown in fig. 3(c). 



2.00299 11.5677 1.69147 1.37704 




0.903849 1.34923 0.373691 0.213635 



U V(°) RetyM) ImtyM) 

FIG. 1: Solutions of the nonlinear problem (3) and the linearized problem (4,5), i.e. the Goldstone modes, at the "best" parameters, p max = 25, 
N p = 1280, Ng = 64, as density plots. Numbers under the density plots are their amplitudes A: white of the plot corresponds to the value A 
and black corresponds to the value —A of the designated field. Upper row: 1st components, lower row: 2d components. 




2.00299 0.0439188 0.72395 0.452911 




0.903849 0.162824 2.70424 2.80215 



U W(°) RefwW) IrnfwW) 

FIG. 2: Same visualization as in fig. 1, for the adjoint linearized problem (13,15), i.e. the response functions. 
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FIG. 3: Radial dependence of the angle-averaged solutions for the spiral wave (a), Goldstone modes (b) and response functions (c). In (c), the 
dependence of Ati;(p max ) = ii(p max ) — d)(25) is shown for comparison, where cD(p max ) is the numerically found spiral angular velocity in 
the disk of given radius p max . 



Sandstede and Scheel [19, 20] have computed exponential 
decay/increase rates of eigenfunctions of periodic wavetrains 
in one spatial dimension. A similar technique should, in prin- 
ciple, also work for the adjoint eigenfunctions. Knowning the 
asymptotic wavelength of the spiral wave, this can be used to 
predict the exponential decay rates of the RFs of spiral waves. 
As can be seen from the results of Wheeler and Barkley [21], 
although such correspondence between ID and 2D calcula- 
tions can be established, the accuracy of decay rate estimates 
for two-dimensional eigenfunctions achieved in this way is in- 
sufficient for a meaningful estimate of the accuracy of those 
eigenfunctions. 



B. Convergence 

We now turn to the main results of our study. Convergence 
of the method has been tested by changing one of the three 
numerical approximation parameters p max , N p and Ng while 
keeping the other two at the fixed values set by the "best ex- 
ample". More specifically, while changing p max , we consider 
two variants: one with fixed Ng, and one with changing Ng 
so that the combination p max A8, which is the size of the out- 
ermost computational cells in the angular direction, remains 
constant. 

Fig. 4 illustrates the results of the study, where the four 
columns correspond to different series of calculations, and the 
three rows correspond to the three different methods of as- 
sessing the accuracy: closeness of the eigenvalues to the theo- 
retical values, distance between "numerical" and "analytical" 
GMs and orthogonality between non-dual RFs and GMs. The 
scales of Ap, A8 and the error estimates are logarithmic, and 
the scales of p max are linear. Here shown is the distance be- 
tween the "numerical" and "analytical" Goldstone modes in 
L2 norm, the distance in Cq norm looks similar. 

A typical feature on many of the curves is a "knee"-shape, 
when the measure of the error decreases as p max grows or A8 
or Ap decrease, but only until a certain point, beyond which 
it reaches a plateau. This behavior is expected and explicable. 
The calculation error is affected by many factors, and if the 
factor varied in a particular series becomes negligible, then the 
error remains at a constant level determined by fixed values of 



other factors. 

The position of the "knees" on the curves indicates that 
the accuracy of the rotational (n = 0) modes would be im- 
proved if A8 were further decreased (there are no knees on 
the curves corresponding to the rotational modes, red online, 
in the fourth, i.e. rightmost column), whereas the limiting pa- 
rameter for the translational (n = 1) modes is Ap (there 
are no knees on the curves corresponding to the tanslational 
modes, blue online, in the third column). The analysis of 
the first two columns is more complicated. The errors esti- 
mates at the maximal p max are similar in both columns as 
they correspond to the same "best" spiral. These limit val- 
ues are achieved, i.e. plateaux are observed, at much smaller 
Pmax values if Ad = const, than if p max A6* = const. This is 
because reduction of p max at fixed AO produces an additional 
improvement of approximation due to the angular discretiza- 
tion. When p max A6* is kept fixed, as in the first column, the 
dependence of the solution on the disk radius is without this 
extra benefit. 

The rates of convergence with respect to parameters can be 
assessed by the slopes of the curves above the knees before 
they plateau. In some cases the data is somewhat irregular, 
primarily at parameters corresponding to lower values of er- 
ror estimates. This is not unexpected and we attribute it to in- 
complete convergence of the iterative procedures (see below). 
On the whole, the slopes can be determined clearly from these 
plots. 

The constant slope in the first (leftmost) and the second 
columns corresponds to the exponential convergence with 
Pmax- The constant slope in the third column corresponds to 
power-law convergence, and the typical slope is 2. This is 
well seen on the curves for translational modes, blue online, 
and not well on the curves for rotational modes, red online, 
which are very small anyway. Slope 2 in the third column 
is to be expected as our discretization is second-order in Ap 
in all cases. The curves in the fourth (rightmost) column are 
convex, which is consistent with the fact that the order of ap- 
proximation is Ng, which varies along the curve as A8 varies, 
since Ng = 2ir/A8, so the slope is bigger for smaller A8. In 
other words, the high order of the Fornberg approximation of 
the 8 derivatives implies the convergence in A8 is faster than 
any fixed power. 



8 




3 10 12 14 16 13 20 22 24 26 3 10 12 14 16 13 20 22 24 26 10" 2 10" 1 10° 0.1 0.2 0.3 



(l) Pm, (RA8 = const) (j) p„„ (A8 = const) (k) Ap (]) A6 

FIG. 4: Convergence in numerical parameters: of deviation of the numerical eigenvalues from theoretical (upper row), of L2 distance between 
numerical and theoretical eigenfunctions (second row) and of orthogonality, i.e. Frobenius norm of the difference of the matrix of scalar 
products of eigenfunctions and adjoint eigenfunctions from the unity matrix (third row), all in logarithmic scales, as dependencies of disk 
radius (first and second columns, linear scale), radius discretization step (third column, logarithmic scale) and polar angle discretization step 
(fourth column, logarithmic scale). In the first column, /5 max is changed while the values of Ap and RAO are kept constant. In the second 
column, pmax is changed while Ap and A6 are kept constant. 
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FIG. 5: (a,b) Effect of the accuracy of the unperturbed spiral wave solution on the convergence: (a) Newton-iteration tolerance 10 . (b) 
Newton iterations until the norm of the residual stopped decreasing, (c) Convergence of the response functions in Ap. 



The irregular shape of some of the curves in fig. 4 at very 
low values of the error estimates is related to the accuracy of 
finding the spiral solution and is ulitmately affected by the 
precision of floating point computations. Note that all cal- 
culations in fig. 4 have been performed with a tolerance of 



10 for Newton iterations of the spiral wave and some of the 
curves fall as low as 10~ 15 i.e. close to machine epsilon. A 
change in the tolerance of the Newton iteration reduces irreg- 
ularities in the curves at low values, as shown in fig. 5(a,b). 

Finally, fig. 5(c) illustrates convergence of numerical RFs 
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W^ 0,1 ' as Ap — > 0, calculated as the L2-distance £0.1 be- 
tween the solutions at a given resultion Ap and the "best" 
solution calculated at the smallest Ap* = 25/1280. As ex- 
plained in the Sec. II B, this comparsion has been restricted 
to the series of calculations with varying Ap, where grids at 
lower resolutions were subgrids of those with higher resolu- 
tions. The graphs of Co distances £' Q x looked similar and are 
not shown here. 



IV. DISCUSSION 

The main result of this paper is a general, robust method for 
obtaining response functions for rigidly rotating spiral waves 
in excitable media with required accuracy. 

We have tested the method on the FitzHugh-Nagumo model 
and we have studied the convergence of spiral wave solutions 
and eigenfunctions, both the Goldstone modes and the re- 
sponse functions, with respect to the numerical approximation 
parameters p max , N p and Ng. The rates of convergence are 
found to agree with the order of approximation and indicate 
the accuracy with which solutions can be found for particular 
numerical parameters. 

The slowest (second-order) convergence is, as expected, in 
the parameter N p . Thus in a typical situation, an improvement 
of accuracy requires, other things being equal, an increase of 
N p , with associated increase in memory and time demands. 
Thus, the most promising avenue of further development of 
the method is via increase of the approximation order of the 
radial derivatives. This is, of course, subject to usual caveat 
that the degree of approximation should be consistent with the 
actual smoothness of the solutions. 

The method used here to solve the eigenvalue problems for 
operators L relies on successive application of transforma- 
tions of L applied to a sequence of vectors, alternating with 
Gram-Schmidt orthogonalization. These are typical ideas, 
also used in [40, 41]. The difference is that in [40, 41], the lin- 
ear transformations were polynomial functions of L whereas 
we use rational functions of L. The polynomial iterations used 
in [40, 41] were in fact equivalent to solving a Cauchy prob- 
lem for equation du/di = Lu by the explicit Euler method. 
Therefore, those methods require a large number of iterations, 
and convergence speed of the iterations depends on the small- 
ness of the absolute difference of the real parts of the eigen- 
values of interest compared to those of other eigenvalues. One 
requires at least O(10 5 ) and typically O(10 6 ) sparse matrix- 
vector multiplications to achieve the desired solutions to the 
eigenvalue problem using such an approach. 

In contrast, with the complex shift and inversion of L used 
in this paper, the convergence speed of the iterations depends 
on the smallness of the distance of the eigenvalues from their 
theoretical values used in the complex shift, compared to the 
distance to other eigenvalues. Hence the number of iterations 
required is very small, typically O(10). More specifically, 
with Krylov subspace dimensionality 3, the number of ma- 
trix multiplications with matrix B of (17) did not exceed 7 
per one eigenpair; with Krylov subspace dimensionality 10, 
this number rose to 10. The price to pay for this accelera- 



tion is the necessity to solve large systems of linear equations. 
However, the key observation is that since the linear system 
is fixed, it needs to be factorized only once, for a given com- 
plex shift, and used for all iterations. Multiplication by ma- 
trix B is achieved with only inexpensive back/forward solves. 
Moreover, due to the way we ordered the unknowns in the 
discretized problem, the sparcity of matrix B does not depend 
on the order of approximation of ^-derivatives. Hence, we are 
able to employ high-order approximations requiring far fewer 
points in the 9 direction for the same accuracy as the second- 
order finite difference discretization used in [40], thereby fur- 
ther improving the efficiency of our method. 

Discounting the factorization step, each iteration, which in- 
volves multiplication by B, is comparable to multiplications 
by L. In practice we find that the factorization itself does 
not require more than the equivalent of four to six actions of 
B. On a MacPro with 3 GHz Intel processor, the factorization 
step takes e.g. about 7.5 sec for the grid N p = 1280, Ng = 64, 
and 0.67 sec for the grid N p = 640, Ng = 32; the compu- 
tation times per B-multiplication were 1.23 sec and 0.17 sec 
respectively. 

The comparison of our present method with [41] is un- 
equivocal: matrix inverses were not used there, and it was 
admitted already in [41] that the resulting accuracy of solu- 
tions was severely limited. While direct accuracy and timing 
comparisons with [40] would be most convincing, that code 
is not publicly available. However, for reasons already noted, 
on any given polar grid, the method we report is more accu- 
rate due to the angular discretization and considerably faster 
in floating-point operations. 

The computed response functions are localized in the vicin- 
ity of the spiral wave tip and exponentially decay with dis- 
tance from it. This localization ensures convergence of the 
convolution integral in (12) in an unbounded domain. 

The eigenvectors of the linearized operator, i.e. Goldstone 
modes and of its adjoint, i.e. the response functions have been 
computed using the same technique, so the qualitatively dif- 
ferent behavior of these solutions at large p is not a numerical 
artefact, as it was not in any way assumed in the numerical 
method. 

Although the method has been used here to compute the re- 
sponse functions in the FitzHugh-Nagumo model, none of the 
details of the method depends on any specifics of the partic- 
ular reaction kinetics and should be widely applicable to the 
computation of response functions of rigidly rotating waves in 
any other model of excitable tissue, as long as its right-hand 
sides are continuously differentiable so the linarized theory is 
applicable. Moreover, the method can also be extended in a 
straightforward way to include additional effects, such as the 
effect of uniform twist along scroll waves with linear filaments 
in three dimensions [17, 40, 48]. 
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